
#Clear R
rm(list=ls())
zero1 <- function(x, minx=NA, maxx=NA){
  res <- NA
  if(is.na(minx)) res <- (x - min(x,na.rm=T))/(max(x,na.rm=T) -min(x,na.rm=T))
  if(!is.na(minx)) res <- (x - minx)/(maxx -minx)
  res
}


#function to install packages if they don't exist
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

# usage
packages <- c("ggplot2", "psych","interplot","stargazer","dplyr","stringr","tidyr","lavaan")
ipak(packages)


#load data
load("Study2_data.Rdta")



#############################
#
#Appendix B3 check hear/see treatment is not predicted by treatment  
#
#############################
systematicdrop<-glm(dumHearSee~t_message+t_partycue+t_gender+t_aggressive,data=nc, family="binomial")
summary(systematicdrop)

stargazer(systematicdrop, title="Check if exclusion of respondents is randomly distributed across conditions", align=TRUE, 
          dep.var.labels=c("Could not hear and/or see treatment (1)"), 
          omit.stat=c("LL","ser","f", "adj.rsq"), 
          covariate.labels = c("Policy Information", "Ideological Cue", "Gender",
                               "Aggressive Speech", "Constant"),
          star.cutoffs=c(0.05), 
          notes = "Logistic Regression models; *p<0.05", 
          notes.append = FALSE,no.space = T)

#############################
#
#subset - exclude those who could not hear/see the treatment
#
#############################
nc <- subset(nc,dumHearSee==0)


#############################
#
# Figure 3
#
#############################

#run models
p1 <-       lm(zero1(dv_agree_partyposition)~comb_NFC_ANES*t_message+t_partycue*comb_NFC_ANES+t_gender*comb_NFC_ANES+t_aggressive*comb_NFC_ANES,subset(nc))
summary(p1)
p2 <-       lm(zero1(dv_agree_partyposition)~comb_NFC_Bullock*t_message+t_partycue*comb_NFC_Bullock+t_gender*comb_NFC_Bullock+t_aggressive*comb_NFC_Bullock,subset(nc))
summary(p2)
p3 <-       lm(zero1(dv_agree_partyposition)~comb_NFC_full*t_message++t_partycue*comb_NFC_full+t_gender*comb_NFC_full+t_aggressive*comb_NFC_full,subset(nc))
summary(p3)
#marginal effects plots
p1inter <- interplot(m = p1, var1 = "t_message", var2 = "comb_NFC_ANES")$data
p2inter <- interplot(m = p2, var1 = "t_message", var2 = "comb_NFC_Bullock")$data
p3inter <- interplot(m = p3, var1 = "t_message", var2 = "comb_NFC_full")$data
p1inter$index <- "ANES 2-item"
p2inter$index <- "Bullock 6-item"
p3inter$index <- "Cacioppo et al. 18-item"

forplot <- data.frame(rbind(p1inter,p2inter,p3inter))
forplot$index <- factor(forplot$index,levels=c("ANES 2-item","Bullock 6-item","Cacioppo et al. 18-item"))
ggplot(forplot, aes(x=fake,y=coef1))+geom_line()+geom_ribbon(aes(ymax=ub,ymin=lb),alpha=.4)+facet_wrap(~index)+theme_bw()+ylab("Marginal Effect of Policy Information")+xlab("Need for Cognition")

###########################################
#
#Figure 4
#
###########################################

#fit 
names(nc)
ncitems <- nc[,81:98]
nclist <- list()
## Create Empty List to stick output into  k = the number of combinations when we have between 2 and 17 items in a set
k <- 0
for(i in 2:17){
  k <- k + ncol(combn(1:18, i))
}

nclist <- list()
# empty matrix to stick alpha scores
alphacoefs <- data.frame(matrix(nrow=k,ncol=4))
h=1
for(j in 2:17){
  # create possible combinations of item numbers of length j
  combos <- combn(1:18, j)
  # empty matrix to stick individual means made of of items from combinations set in combos
  
  two <- matrix(nrow=nrow(nc),ncol=ncol(combos))
  for(i in 1:ncol(combos)){
    two[,i] <-   rowMeans(ncitems[,combos[,i]])
    alphacoefs[h,] <-   psych::alpha(ncitems[,combos[,i]],check.keys = TRUE)$total[1:4]
    h=  h+1
  }
  ## m
  
  lmcoefs <- matrix(nrow=ncol(combos),ncol=4)
  for(i in 1:ncol(two)){
    lmcoefs[i]   <-     lm(zero1(nc$dv_agree_partyposition)~zero1(two[,i])*nc$t_message+nc$t_partycue)$coefficients[5]
  }
  print(j)
  nclist[[j]] <- lmcoefs
}

b <- 0
items <- list()
for(j in 2:17){
  b=b+1 
  combos <- combn(1:18, j)
  for(i in 1:ncol(combos)){
    items[[b]] <- combos
  }}

itemrow <- lapply(items,nrow)
itemcol <- lapply(items,ncol)

collapsepastecolumns <- function(x)apply(x,2,paste,collapse=",")
rr <- lapply(collapsepastecolumns,X = items)

numbero <-   unlist(itemcol)
emptylist <- list()  
for(i in 2:17){
  emptylist[[i-1]] <- rep(i,numbero[[i-1]])
}




coefficients <- do.call("rbind",nclist)
itemstogether <- unlist(rr)
numall <- unlist(emptylist)


library(ggplot2)
both <- data.frame(coefficients,alphacoefs,numall,itemstogether)

eighteen <- rowMeans(ncitems,na.rm=T)
bottomrow <- summary(lm(zero1(nc$dv_agree_partyposition)~zero1(eighteen)*nc$t_message+nc$t_partycue))$coefficients[5,]
names(bottomrow) <- names(both)[1:4]

oo <- rbind(both,NA)
oo[nrow(oo),1:4] <- bottomrow
oo[nrow(oo),]$numall <- 18

for(i in 1:18){
  bottomrow <- summary(lm(zero1(nc$dv_agree_partyposition)~zero1(ncitems[,i])*nc$t_message+nc$t_partycue))$coefficients[5,]
  oo <- rbind(oo,NA)
  oo[nrow(oo),1:4] <- bottomrow
  oo[nrow(oo),]$numall <- 1
}


head(alphacoefs)



save(oo,file="study2randomcoefsalpha.RData")

ggplot(oo,aes(y=X1.1,x=as.factor(numall)))+geom_boxplot()+theme_bw()+xlab("Number of items")+ylab("Coefficient")

ggplot(subset(oo,numall!=1 & numall!=18),aes(y=X1.1,x=X1))+geom_point()+theme_bw()+xlab("Cronbach's Alpha")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_wrap(~numall)

## split up which items went into which coefficient
load("study2randomcoefsalpha.RData")
oo %>% separate(itemstogether, into = paste("V", 1:18, sep = "_")) -> both

CF18<-'NFC_full18 = ~ nfc1+ nfc2+nfc3+nfc3+nfc4+nfc5+nfc6+nfc7+nfc8+nfc8+nfc9+nfc10+nfc11+nfc12+nfc13+nfc14+nfc15+nfc16+nfc17+nfc18
# fix variance of speed factor
NFC_full18 ~~ 1*NFC_full18'

oo %>% separate(itemstogether, into = paste("V", 1:18, sep = "_")) -> both



itemscores <- inspect(fit,what="std")$lambda
itemsforfactor <- dplyr::select(both,starts_with("V_"))

for(i in 1:18){
  itemsforfactor[itemsforfactor==i] <- itemscores[i]
}

itemsforfactor <- mapply(as.numeric,itemsforfactor)
factormeans <- rowMeans(itemsforfactor,na.rm=T)
both$numall_1 <- factor(paste(both$numall,"Items"))
both$numall_1 <-  factor(both$numall_1,levels=levels(both$numall_1)[c(1,11:18,2:9)])
all <- data.frame(both,factormeans=factormeans)


ggplot(subset(all,as.numeric(numall)!=1  & as.numeric(numall)!=18),aes(y=X1.1,x=factormeans))+geom_point()+theme_bw()+xlab("Mean Factor Score of Items")+geom_smooth(method="lm")+ylab("b Coefficient")+facet_wrap(~numall)

